毕业季，把毕设某Goldstein算法和FFT最小二乘算法实现后就扔了，借了一大堆逆向汇编调试的书。

我想让大家看看我是如何瞎折腾最终浪费生命一事无成的，请大家从这个反例中吸取经验教训，珍爱生命，远离折腾。

专心赚钱和把妹，别搞些没用的。

中国谚语

Gray Hat Python

跳坑从一本Python灰帽子开始，因为对Python的熟悉和兴趣加上某种机缘巧合，从图书馆借来看了一阵。

首先是开发环境的搭建，据作者说搭建开发环境是个痛苦而无聊的过程。对新手很棘手、很消磨斗志和耐心。但对我这样已经在Linux下被虐过多少遍的人完全没感觉了，即使在windows上也轻车熟路装好Python。唯一不爽的地方在于，作者说windows下调试乐趣横生所以不考虑linux，我只好装了个虚拟机xp玩。当然，作为有个性的Linuxer，我的做法是在linux宿主机编辑文件共享给windows虚拟机用，真尼玛蛋疼。

接着就是些基本理论，通过这本书我第一次发现原来Python和C集成的这么好，于是把Python C API翻译了一遍……后来在《Unix设计艺术》中看作者说Python最开始设计就是为了跟C集成……

总之，Python强大的C外挂= =。但总有些不舒服的地方，尼玛Windows和Unix函数参数进栈顺序……

再往后是调试器设计的基本知识，卧槽啊，第一次了解CPU，那个看得一愣一愣的，几次差点睡着。最后看了好几遍又回头了好几遍又看了好几本其它书讲才总算知道有几个寄存器且有什么用栈是什么软硬件断点实现，现在什么平坦内存模型内存断点堆是什么还是搞不清，传说中的Intel开发者手册更是没勇气去看。

接着就开始在虚拟机xp上开始照着自己写调试器了，好在作者对各种Win32 API讲得清清楚楚明明白白，又在网上提供源码和各种结构定义。不然您看着满目错误的中文版就跪吧。照着网上提供的代码试了不知道多久终于把这个向pygdb致敬的调试器写出来了。

不过Win32 API完全给我看跪了。在byr上看有人说Linux的系统调用比windows好懂[^1]。最开始还想将这个具简单的调试器移植到linux下看看，然而，现在还什么都没看。

之后，慢慢看吧，还没动手，用pygdb部署软钩子，进程快照和回滚，处理非法内存操作貌似很厉害的样子……

再往后看到用Immudbg针对特定shellcode寻找漏洞，我觉得得回头找点其它书补补相应知识了，很难看懂作者在说什么了……于是，IDA pro

IDA Pro

前面讲汇编和IA-32处理器，栈及堆的基本原理可执行文件的区段指令集啥的，看得云里雾里，不知道是因为看雪翻译的问题还是我弱爆了……

看到PE格式和ELF格式文件格式细节完全跪了没，最终也没理解什么GOT/PLT，动态链接成员结构啥的链接加载啥的……然后就跳过了。于是后来看到《Debug Hacks》听作者讲GOT/PLT时又跪了……

接着看了看第一个例子，跟踪执行流。虽然下了一个IDA Pro以供调教，然而终于也没开虚拟机……我说，我还是回头看看汇编在来看例子吧……于是这本书又暂时按下。

Debug Hacks

一本日本人写的书，风格清新脱俗，平易近人。看讲gdb那货直接把我迷住了，哇草啊，gdb这么好玩啊，哇草啊，比我写的那啥调试器强大多了。

看到byr挂站排行榜脑子一抽写了个byr挂站脚本扔到Luoyun云平台上，结果web-ssh关掉了看不到标准输出了。然后竟然在Stack Overflow上看到有人用gdb来重定向程序输出，尼玛啊，这瞬间就感觉调试器无敌了。

结果就是，在看到这本主要讲linux平台调试的书时顿时爱不释手，特别是前几章讲gdb的，再加上几个例子太好玩了。

最sb的事接着来了，参照书上写了个对NULL写入的Linux内核模块接着竟然在宿主机上加载了……卧槽，久违的kernel panic。好吧，正好复习下reisub大法，可是它没重启动，最后还是长按关机键……

后面讲crash倒兴致了然了，虽然在gentoo上顺手装上，却没啥动力让系统吐核……虚拟机又懒得开。

顺口吐槽，原来pdb这么像gdb啊……

不过还是好多不懂，linux下线程共享资源，看得完全跪掉了，现在也没搞明白怎么回事。

gdb或objdump反汇编出来的汇编代码语法是AT&T的，和Intel差别还挺大，之前讲windows的IDA pro和Gray Hat Python都是Intel语法……我喜欢AT & T

于是抛弃了几本讲MSMA和NSMA的书选择了《Professional Assembly Language》。

Professional Assembly Language

一本完全遵循AT&T语法，在linux平台下，实用gcc/gdb/ld/as/gas/grof等工具的书啊。

第一个例子无压力通过，虽然不了解语法云里雾里，但之前硬件基础讲解还是很棒的，就是只是32位的，我的系统是64兼容32位，有时候要加些参数。相比而言，在这个问题上只有《Debug Hacks》做的好，64位和32位都讲得挺清楚。

打算接着熟悉下指令

再接着打算看看Linux系统调用。回头学习Linux去= =

以上只是打算……

K & R

为什么乱入了一本书，带着浓重Unix风格的C语言圣经。为了看懂那些书中的例子在调试些什么C程序，只好翻K & R去了。

K & R也许真不适合入门，但我还是想给初学者推荐(莫喷……)。

听说K & R不符合C 99规范，暂不管它，Fortran这规范都2008了也没听说哪个编译器实现了，还有好多东西还是Fortran77。gcc不报Warning就好。

Coursera

最不能理解的是，这个时候还跑去看High Performance Scientific Computing和Machine Learning。开着八线程学Fortran和Octave(matlab)。

Summary

大概就是这样，被我浪费掉的，这么多人都忙着旅行恋爱的享受最后青春的，我选择的，大学最后光阴。

大家跟我一起唱：

珍爱生命，远离计算机
珍愛生命，遠離計算機
Cherish life, away from the computer
大切な命のため、コンピュータから離れてください
Chérir la vie, loin de l'ordinateur
Берегите жизнь, вдали от компьютера
Cherish Leben, weg vom Computer
컴퓨터에서 멀리, 생명을 소중히
Apreciar la vida, lejos de la computadora
++++++++[>+>++>+++>++++>+++++>++++++>+++++++>++++++++>+++++++++>++++++++++>+++++++++++>++++++++++++>+++++++++++++>++++++++++++++>+++++++++++++++>++++++++++++++++<<<<<<<<<<<<<<<<-]>>>>>>>>+++.---<<<<<<<<>>>>>>>>>>>>>.<<<<<<<<<<<<<>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<>>>>>>>>>>>>>>++.--<<<<<<<<<<<<<<>>>>>>>>>>>>>+.-<<<<<<<<<<<<<>>>>>>>>>>>>>>+++.---<<<<<<<<<<<<<<>>>>>>>>>>>>>.<<<<<<<<<<<<<>>>>.<<<<>>>>>>>>>>>>>>----.++++<<<<<<<<<<<<<<>>>>>>>>>>>>>+.-<<<<<<<<<<<<<>>>>>>>>>>>>>--.++<<<<<<<<<<<<<>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<>>>>>>----.++++<<<<<<>>>>.<<<<>>>>>>>>>>>>+.-<<<<<<<<<<<<>>>>>>>>>>>>>>>-.+<<<<<<<<<<<<<<<>>>>>>>>>>>>+.-<<<<<<<<<<<<>>>>>>>>>>>>>>>+.-<<<<<<<<<<<<<<<>>>>.<<<<>>>>>>>>>>>>>--.++<<<<<<<<<<<<<>>>>>>>>>>>>>>++.--<<<<<<<<<<<<<<>>>>>>>>>>>>>>-.+<<<<<<<<<<<<<<>>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<<>>>>.<<<<>>>>>>>>>>>>>>>----.++++<<<<<<<<<<<<<<<>>>>>>>>>>>>>.<<<<<<<<<<<<<>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<>>>>.<<<<>>>>>>>>>>>>+++.---<<<<<<<<<<<<>>>>>>>>>>>>>>-.+<<<<<<<<<<<<<<>>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<<>>>>>>>>>>>>>>.<<<<<<<<<<<<<<>>>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>----.++++<<<<<<<<<<<<<<<>>>>>>>>>>>>>---.+++<<<<<<<<<<<<<.

人生得意须尽欢：）

花堪折时只须折：）

补记

Fri 24 May 2013 07:03:04 PM CST

谣传要交论文，基本上是两天给弄完了，好在之前真是认真了解这个课题了，吹起牛来一气呵成。

快毕业了一堆事，又是转党组织关系的又是填毕业登记表，又是毕业照体检。

班里没有毕业旅行感觉太好了，毫无兴致。

毕业把妹这种事……我还是不抱希望了……

整日整日地沉迷在图书馆和计算机前，偶尔的偶尔一个人出去转转。

哈，东边小操场没什么人，太喜欢在那边跑步了，江南的夜色依旧很美。记得有次跟着一个妹子跑了半个小时，第二天直接虚掉了= =

这几天就看了Professional Assembly Language一本书，断断续续看到浮点计算。把书上的代码跟着打了一遍，不得不说这本书真棒。例子都是神来之笔。

开始看《Hacking: The Art of Exploitation》, C语言那部分的例子真是让人觉得相见恨晚，由于看过以上几本书对汇编和堆栈比较熟悉发现看着真是享受。准备把书中所有的代码打一遍，太好玩了。

借了本《别闹了，费曼先生》，看得五体投地，这才是真黑客，真物理学家和数学家。

借了本《数学大师：从芝诺到庞加莱》，不得不说多好的一本书译得真有些不好懂……

这两本书我想单独写读后感。

因为毕业季一堆事，Coursera上Machine Learning这周到现在基本没看，High performance Scientific Computing直接被落下了，上周的视频还没看完，再加上进入Fortran和makefile的课题了，总有些力不从心。

最后吐槽下Girl-atlas，新上传的图像不仅打上水印而且只能查看缩略图了，而且下载链接也被移除了。很不爽，不能做壁纸的小图完全没有兴致再用脚本抓了。

结果愤而转向其它网站开始抓图……比如Motu( Maybe 18X,慎点 ), 再也没有girl-atlas老少咸宜高清无码的品味了。无论是图像更新及时度，图像质量，网站质量，girl-atlas做的都是最棒的.

FootNotes

[^1]: Deepin 开发者lovesnow(小邪兽)在贴吧贴过如何枚举linux下进程，好厉害好厉害，以前还见过他用python来hack植物大战僵尸= =，膜拜膜拜膜拜。

转眼又是一年的五月中旬，我都不能记起之前的三年五月是什么样，我干了什么，想些什么，爱过什么，恨过什么。

只是这个五月，笼罩在无所事事的清闲和对未来捉摸不定的焦虑气氛中。

毕业论文

导师说，五月末就交终稿了。

然而一个字也不想写。

下雨的某天和室友去教室写了点论文，把绪论写完又写了一堆原理。相比一群复制粘帖的同学我竟然一直能自己有的写……

只是Goldstein算法实践了之后再也不想做什么了。

Coursera

还有就是跟Coursera上课程，总有人喜欢做些 没用 的事。我觉得真该找些家教兼职出去丰富下社会经验换点外快，然而，总有些残想，唉。

认认真真地开始跟Coursera上的machine learning课程和High Performance Scientific Computing. 后来看到Introduction to data science和一门类似营养与健康的课。

最后，Introduction to data science作业难度实在够呛，营养健康感觉太欧美化了。最终只跟了前两门课。结果High Performance Scientific Computing课程在讲到使用Amazon Aws云服务时我跪了，尼玛屌丝没有信用卡……

Machine Learning在Coursera上比较简单，略过了让人抓狂的数学推导而集中精力讲解怎么实际实现算法和解决问题，主要用matlab或octave。感觉相当顺手唉。

High Performance Scientific Computing则讲了些我很熟悉的IPython + numpy +scipy科学计算工具集，还讲了点linux和Python，一点git的用法。总之，太熟悉了实在没什么压力。

Introduction to Data Science有的video没字幕，那老师的吐字我直接就跪了，后来发现有字幕的也不怎么看得懂……后来做了Problem012之后觉得作业压力大，几门课又是毕业论文还有其它想法，还是放弃这门课了。

最后门营养与健康的课。没想到上来就被各种计量单位和机构名称给搞跪了。

爬虫

High Performance Scientific Computing课程的视频不允许下载，这很奇怪，也不合Coursera世界上每个人都能自由获取教育的理念。所以，我找了个冠冕堂皇的理由写了个脚本把用来下相应的MP4视频和字幕。

看到byr挂站排行榜写了个挂站脚本，扔到了Luoyun上。

重新装上了selenium。有空看看Grobot。

其它

开着虚拟机写word论文，写着写着闻到糊味了……

后来zy就把笔记本散热器借我了。想论文写完自己拆机。

在图书馆看了看《漫画SQL》，果然只有漫画看得进去。

发现Gdb真是强大啊，今天竟然用来重定向其它pts程序的标准输出了，这些日子又零零散散看些汇编和逆向的基本知识，感觉又误入歧途了。

ps:难得还打了两次网球几次乒乓球。

pps: 图书馆妹子真多。

ppps：可惜不关我什么事，万花丛中过，片叶未沾身。

pppps：毕业旅行……毫无兴致呢

晚上又一个人独享了JNU摄人心魄的绝美夜色。从操场向东远望城市灯火映红天空，一轮初生的圆月，像一枚猩红的宝石，静悄悄凝视着大地。好像用画笔把那份惊艳记录下来。

这样的春夜将不多， JNU剩下的每个日子都将弥足珍贵。

回到学校这段日子，也许将是‘空前觉后’的安逸日子。然而这安逸平静的日子却弥漫着难抑的忧伤，朦朦胧胧，好像江南氤氲的湿气一样，即使在明媚祥和的春日里，也渗透得无处不再。

十六号回学校，今天二十七，除了些零零碎碎的事(比如面基。。)这十一天我开了几个线程：

• 回来在导师面前晃了下，表示我回学校了，又问了问发现毕设难度又被降了一个等级。又把那之前在家写的不堪入目的添加branch cuts的七层for循环重新调试好了，结果今天发现去包裹算法有问题，竟然会出现2 pi整周期偏差，好吧，在老师面前还吹嘘简单……^1
• 当我看到girl-atlas这个有魅力的网站时，第一反应是把喜欢的写真集下载下来。结果这站用的竟然是什么城通网盘，翻越教育网的代理环境下根本无法下载。然后……竟然习惯性的打开page source细细查看，然后很熟练的写起来脚本抓图……后来，就诞生了这货girl-atlas.py，用来断断续续下了2.36G的写真图片。
• 开始看图书馆里借来的《python灰帽子》,饶有兴趣地在虚拟winxp下动手用Python实现一个简单的调试器。我就不说这本书错误多多了，算是明白豆瓣评分为什么这么低，名副其实……一行一行地敲代码并且调试N遍使之通过测试真不容易。

大概就这，看到了有些东西，有空想折腾：

• 让我把零零碎碎地把js语法赶紧看完
• verletjs - 物理仿真js库
• flask， 对girl-atlas的图片有些想法 XD

FootNotes

翻译自scipy lecture notes: Interfacing with C

作者：Valentin Haenel

这个章节包含许多在python代码中支持c/c++本机代码的许多不同方法， 通常这个过程叫作包裹(wrapping)。本章的目的是让您大致知道有那些技术和它们分别的优缺点是什么，于是您能够为您自己的特定需要选择何时的技术。在任何情况下，一旦您开始包裹，您几乎一定将想要查阅您所选技术各自的文档。

目录

• toc
  {: toc}

简介

本章节包含以下技术：

• Python-C-Api
• Ctypes
• SWIG(简单包裹和接口生成)
• Cython

这四种方法大概是最著名的，其中Cython可能是最高级且应该优先使用的。如果您想从其它角度理解包裹问题，其它方法也很重要。已经说过，虽然还有其它方法，但是理解以上基本方法，您将能评估您自己的选择看是否符合自己的需要。

以下标准在评估一项技术时也许有用：

• 是否需要额外的库？
• 代码是否是自动生成的？
• 需要编译吗？
• 和Numpy数组进行交互方便吗？
• 支持C++吗？

首先，您应该考虑你的用例。当用本机代码接口时，通常有两个用例：

• 存在需要充分利用的C/C++代码，或者那些代码已经存在，或者那些代码更快。
• Python代码慢爆了，将内循环交给本机代码处理

每个技术通过包裹math.h中的cos函数实现。尽管这是微不足道的例子，它将很好的展示基本的包裹问题。因为每个技术也包括某种形式的Numpy支持，这也通过使用一个余弦函数被在某种数组上计算的例子来展示。

最后但重要的是两个小警告：

• 所有这些技术都可能造成Python解释器崩溃(段错误)，这(通常)是C代码的Bug。
• 所有例子在linux上完成，也应该能在其它操作系统上实现
• 大多数例子中你需要一个C编译器

Python-C-Api

Python-C-API是标准Python解释器(就是所谓的CPython)的支柱。使用这个API可以用C或C++语言编写Python扩展。显然这些扩展模块可以凭借语言兼容性，调用任何C或C++写成的函数。

当使用Python-C-API时，人们通常写许多样板代码，先解析传递给函数的参数，然后构建并返回类型。

优点

• 无需额外的库
• 许多低级的控制
• 完全可以用C++

劣势

• 可能需要大量工作
• 代码中的大量开销
• 必须编译
• 高额的维护代价
• 当跨Python版本时若C-Api变化没有后向兼容性

注意：以下Python-C-Api示例主要为了展示需要。因为大多其它技术实际上依赖这个，所以最好对它如何工作有个高层次的了解。在99%的用例中你最好使用其它技术。

示例

以下C扩展模块，让标准数学库中的cos函数在Python中可用：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

/*  Example of wrapping cos function from math.h with the Python-C-API. */

#include <Python.h>
#include <math.h>

/*  wrapped cosine function */
static PyObject* cos_func(PyObject* self, PyObject* args)
{
    double value;
    double answer;

    /*  parse the input, from python float to c double */
    if (!PyArg_ParseTuple(args, "d", &value))
        return NULL;
    /* if the above function returns -1, an appropriate Python exception will
     * have been set, and the function simply returns NULL
     */

    /* call cos from libm */
    answer = cos(value);

    /*  construct the output from cos, from c double to python float */
    return Py_BuildValue("f", answer);
}

/*  define functions in module */
static PyMethodDef CosMethods[] =
{
     {"cos_func", cos_func, METH_VARARGS, "evaluate the cosine"},
     {NULL, NULL, 0, NULL}
};

/* module initialization */
PyMODINIT_FUNC

initcos_module(void)
{
     (void) Py_InitModule("cos_module", CosMethods);
}

如您所见，所有对参数处理、返回类型和模块初始化都相当样板化。然而有些被摊销了，当扩展增长时，样板需要每个函数保留。

标准python构建系统distutils支持从setup.py编译C扩展，这相当方便。

1
2
3
4
5
6
7

from distutils.core import setup, Extension

# define the extension module
cos_module = Extension('cos_module', sources=['cos_module.c'])

# run the setup
setup(ext_modules=[cos_module])

这能被编译：

1
2
3
4
5
6
7
8
9

 ~/Work/scipy-lecture-notes/interfacing-with-c ⮀ ls
cos_module.c  setup.py
 ~/Work/scipy-lecture-notes/interfacing-with-c ⮀ python setup.py build_ext --inplace
running build_ext
building 'cos_module' extension
x86_64-pc-linux-gnu-gcc -pthread -fPIC -I/usr/include/python2.7 -c cos_module.c -o build/temp.linux-x86_64-2.7/cos_module.o
x86_64-pc-linux-gnu-gcc -pthread -shared build/temp.linux-x86_64-2.7/cos_module.o -L/usr/lib64 -lpython2.7 -o /home/lyy/Work/scipy-lecture-notes/interfacing-with-c/cos_module.so
 ~/Work/scipy-lecture-notes/interfacing-with-c ⮀ ls
build  cos_module.c  cos_module.so  setup.py

• build_ext是用来构建扩展模块的
• --inplace将编译好的扩展模块输出到当前文件夹

文件cos_module.so包含编译的扩展，我们能将它加载到IPython解释器中：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

In [1]: import cos_module

In [2]: cos_module?
Type:       module
String Form:<module 'cos_module' from 'cos_module.so'>
File:       /home/lyy/Work/scipy-lecture-notes/interfacing-with-c/cos_module.so
Docstring:  <no docstring>

In [3]: dir(cos_module)
Out[3]: ['__doc__', '__file__', '__name__', '__package__', 'cos_func']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0

现在让我们看看它有多健壮：

1
2
3
4
5
6
7

In [7]: cos_module.cos_func('foo')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-11bee483665d> in <module>()
----> 1 cos_module.cos_func('foo')

TypeError: a float is required

Numpy支持

类似于Python-C-API,Numpu自身作为C扩展实现也有Numpy-C-API。这个API可以在写自定的C扩展时，被用来从C创建和操作Numpy数组。参见高级Numpy

以下例子展示了如何将Numpy数组作为参数传递给函数，如何使用(老的)Numpy-C-API遍历整个Numpy数组。它仅仅将数组作为参数，运用来自math.h中的余弦函数，并且返回一个新的结果数组。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

/*  Example of wrapping the cos function from math.h using the Numpy-C-API. */

#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>

/*  wrapped cosine function */
static PyObject* cos_func_np(PyObject* self, PyObject* args)
{

    PyArrayObject *in_array;
    PyObject      *out_array;
    PyArrayIterObject *in_iter;
    PyArrayIterObject *out_iter;

    /*  parse single numpy array argument */
    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array))
        return NULL;

    /*  construct the output array, like the input array */
    out_array = PyArray_NewLikeArray(in_array, NPY_ANYORDER, NULL, 0);
    if (out_array == NULL)
        return NULL;

    /*  create the iterators */
    /* TODO: this iterator API is deprecated since 1.6
     *       replace in favour of the new NpyIter API */
    in_iter  = (PyArrayIterObject *)PyArray_IterNew((PyObject*)in_array);
    out_iter = (PyArrayIterObject *)PyArray_IterNew(out_array);
    if (in_iter == NULL || out_iter == NULL)
        goto fail;

    /*  iterate over the arrays */
    while (in_iter->index < in_iter->size
            && out_iter->index < out_iter->size) {
        /* get the datapointers */
        double * in_dataptr = (double *)in_iter->dataptr;
        double * out_dataptr = (double *)out_iter->dataptr;
        /* cosine of input into output */
        *out_dataptr = cos(*in_dataptr);
        /* update the iterator */
        PyArray_ITER_NEXT(in_iter);
        PyArray_ITER_NEXT(out_iter);
    }

    /*  clean up and return the result */
    Py_DECREF(in_iter);
    Py_DECREF(out_iter);
    Py_INCREF(out_array);
    return out_array;

    /*  in case bad things happen */
    fail:
        Py_XDECREF(out_array);
        Py_XDECREF(in_iter);
        Py_XDECREF(out_iter);
        return NULL;
}

/*  define functions in module */
static PyMethodDef CosMethods[] =
{
     {"cos_func_np", cos_func_np, METH_VARARGS,
         "evaluate the cosine on a numpy array"},
     {NULL, NULL, 0, NULL}
};

/* module initialization */
PyMODINIT_FUNC

initcos_module_np(void)
{
     (void) Py_InitModule("cos_module_np", CosMethods);
     /* IMPORTANT: this must be called */
     import_array();
}

我们仍可使用distutils编译这个。然而，我们必须通过使用numpy.get_include()保证包含了Numpy头文件。

1
2
3
4
5
6
7
8
9

from distutils.core import setup, Extension
import numpy

# define the extension module
cos_module_np = Extension('cos_module_np', sources=['cos_module_np.c'],
                          include_dirs=[numpy.get_include()])

# run the setup
setup(ext_modules=[cos_module_np])

为确信它确实能用我们做以下测试脚本：

1
2
3
4
5
6
7
8

import cos_module_np
import numpy as np
import pylab

x = np.arange(0, 2 * np.pi, 0.1)
y = cos_module_np.cos_func_np(x)
pylab.plot(x, y)
pylab.show()

结果将如下图

Ctypes

Ctypes是一个Python的外部函数库。它提供了兼容C的数据类型。并且允许调用DLL或共享库中的函数。它能够被用来将这些库用纯Python包裹。

优势

• Python标准库的一部分
• 不必编译
• 完全用Python包裹代码

劣势

• 需要将代码包裹作为共享库获得(粗略地说就是Windows中的*.dll、Linux下的*.so和Mac OSX的*.dylib)
• 对C++支持不好

示例

如上所述，包裹的代码是纯Python的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

""" Example of wrapping cos function from math.h using ctypes. """

import ctypes
from ctypes.util import find_library

# find and load the library
libm = ctypes.cdll.LoadLibrary(find_library('m'))
# set the argument type
libm.cos.argtypes = [ctypes.c_double]
# set the return type
libm.cos.restype = ctypes.c_double


def cos_func(arg):
    ''' Wrapper for cos from math.h '''
    return libm.cos(arg)

• 寻找和加载库可能依赖于不同的操作系统，检查文档获取细节。
• 这稍微有点虚幻，因为在系统上已经存在编译好的数学库。如果你将包裹一个在自己的库，你将不得不先编译它。这也许需要也许不需要额外的工作。

我们现在如前述那样使用它：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

In [1]: import cos_module

In [2]: cos_module?
Type:       module
String Form:<module 'cos_module' from 'cos_module.py'>
File:       /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/ctypes/cos_module.py
Docstring:  <no docstring>

In [3]: dir(cos_module)
Out[3]:
['__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 'cos_func',
 'ctypes',
 'find_library',
 'libm']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0

正如之前的例子，这个代码稍微健壮一些。尽管错误信息不怎么有用，因它并没告诉我们应该是什么类型。

1
2
3
4
5
6
7
8
9
10
11
12

In [7]: cos_module.cos_func('foo')
---------------------------------------------------------------------------
ArgumentError                             Traceback (most recent call last)
<ipython-input-7-11bee483665d> in <module>()
----> 1 cos_module.cos_func('foo')

/home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/ctypes/cos_module.py in cos_func(arg)
     12 def cos_func(arg):
     13     ''' Wrapper for cos from math.h '''
---> 14     return libm.cos(arg)

ArgumentError: argument 1: <type 'exceptions.TypeError'>: wrong type

Numpy支持

Numpy包含一些对ctypes接口的支持。特别是有导出Numpy数组作为ctypes数据类型的某一属性的支持，并且有将C数组和Numpy数组互相转化的函数。

更多信息参考Numpy Cookbook中相应章节和numpy.ndarray.ctypes和numpy.ctypeslib的API文档。

在以下例子中，让我们考虑一个库中的C函数，这个函数接受一个数组作为输入并输出一个数组，计算输入数组的正弦值并将结果存储在输出数组中。

这个库包含以下头文件(尽管就这个例子不严格需要，为完整性需要我们列出它)：

1

void cos_doubles(double * in_array, double * out_array, int size);

这个实现在C源码中如下：

1
2
3
4
5
6
7
8
9
10

#include <math.h>

/*  Compute the cosine of each element in in_array, storing the result in
 *  out_array. */
void cos_doubles(double * in_array, double * out_array, int size){
    int i;
    for(i=0;i<size;i++){
        out_array[i] = cos(in_array[i]);
    }
}

因为这个库是纯C的，我们不能使用distutils来编译它。必须同时使用make和gcc:

1
2
3
4
5
6
7
8
9
10

m.PHONY : clean

libcos_doubles.so : cos_doubles.o
	gcc -shared -Wl,-soname,libcos_doubles.so -o libcos_doubles.so cos_doubles.o

cos_doubles.o : cos_doubles.c
	gcc -c -fPIC cos_doubles.c -o cos_doubles.o

clean :
	-rm -vf libcos_doubles.so cos_doubles.o cos_doubles.pyc

我们接着可以将之编译到共享库libcos_double.so中(linux下):

1
2
3
4
5
6
7
8

$ ls
cos_doubles.c  cos_doubles.h  cos_doubles.py  makefile  test_cos_doubles.py
$ make
gcc -c -fPIC cos_doubles.c -o cos_doubles.o
gcc -shared -Wl,-soname,libcos_doubles.so -o libcos_doubles.so cos_doubles.o
$ ls
cos_doubles.c  cos_doubles.o   libcos_doubles.so*  test_cos_doubles.py
cos_doubles.h  cos_doubles.py  makefile

接着我们能继续通过ctypes库对(某些类型)Numpy数组的直接支持包裹这个库了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

""" Example of wrapping a C library function that accepts a C double array as
    input using the numpy.ctypeslib. """

import numpy as np
import numpy.ctypeslib as npct
from ctypes import c_int

# input type for the cos_doubles function
# must be a double array, with single dimension that is contiguous
array_1d_double = npct.ndpointer(dtype=np.double, ndim=1, flags='CONTIGUOUS')

# load the library, using numpy mechanisms
libcd = npct.load_library("libcos_doubles", ".")

# setup the return typs and argument types
libcd.cos_doubles.restype = None
libcd.cos_doubles.argtypes = [array_1d_double, array_1d_double, c_int]


def cos_doubles_func(in_array, out_array):
    return libcd.cos_doubles(in_array, out_array, len(in_array))

• 注意连续单维Numpy数组的固有限制，因为C函数要求这种缓冲区。[^1]
• 注意输出数组必须预先分配，例如通过numpy.zeros()，这个函数将写进它的缓冲区。
• 尽管cos_doubles函数的原始参数是ARRAY, ARRAY, int，最终的cos_doubles_func仅仅接受两个Numpy数组作为参数。

像之前一样，我们相信它能够工作：

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pylab
import cos_doubles

x = np.arange(0, 2 * np.pi, 0.1)
y = np.empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab.plot(x, y)
pylab.show()

SWIG

SWIG, 简化包裹接口生成器，是一个将不同高级编程语言包括Python链接到用C和C++写的程序上的软件开发工具。SWIG重要的功能是，它能自动为你生成包裹代码。这就开发时间来说是个优势，也可能是个负担。生成文件趋于巨大，读起来不友好，包裹过程的结果就是多个间接层，可能有点难以理解。

注意：自动生成的C代码使用Python-C-Api。

优势

• 可以自动包裹给定头文件的整个库
• 对C++工作很好

劣势

• 自动生成巨大的文件
• 若出错难以调试
• 陡峭的学习曲线

示例

让我们假设我们的cos函数位于用C写成的cos_module中，源代码文件为cos_module.c。

1
2
3
4
5

#include <math.h>

double cos_func(double arg){
    return cos(arg);
}

头文件为cos_module.h：

1

double cos_func(double arg);

我们的任务是将cos_func暴露给Python。为了用SWIG实现这个，我们必须写一个包含SWIG指令的接口文件。

1
2
3
4
5
6
7
8
9
10
11

/*  Example of wrapping cos function from math.h using SWIG. */

%module cos_module
%{
    /* the resulting C file should be built as a python extension */
    #define SWIG_FILE_WITH_INIT
    /*  Includes the header in the wrapper code */
    #include "cos_module.h"
%}
/*  Parse the header file to generate wrappers */
%include "cos_module.h"

如您所见，需要太多代码了。在这个简单的例子中在接口文件中仅仅包含头文件就足够将函数暴露给Python。然而，SWIG允许更细粒度地包含/排除头文件中的函数，查看文档获取更多细节。

产生编译的包裹代码是一个两个阶段的过程：

1. 对接口文件运行swig生成文件cos_module_wrap.c，这是用来自动生成Python的C扩展的源代码文件。cos_module.py是自动生成的纯Python模块。
2. 编译cos_module_wrap.c为_cos_module.so。幸运的是，distutils知道如何处理SWIG接口文件，所以我们的setup.py很简单：

1
2
3
4

from distutils.core import setup, Extension

setup(ext_modules=[Extension("_cos_module",
      sources=["cos_module.c", "cos_module.i"])])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

$ cd advanced/interfacing_with_c/swig

$ ls
cos_module.c  cos_module.h  cos_module.i  setup.py

$ python setup.py build_ext --inplace
running build_ext
building '_cos_module' extension
swigging cos_module.i to cos_module_wrap.c
swig -python -o cos_module_wrap.c cos_module.i
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/include/python2.7 -c cos_module.c -o build/temp.linux-x86_64-2.7/cos_module.o
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/include/python2.7 -c cos_module_wrap.c -o build/temp.linux-x86_64-2.7/cos_module_wrap.o
gcc -pthread -shared build/temp.linux-x86_64-2.7/cos_module.o build/temp.linux-x86_64-2.7/cos_module_wrap.o -L/home/esc/anaconda/lib -lpython2.7 -o /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/swig/_cos_module.so

$ ls
build/  cos_module.c  cos_module.h  cos_module.i  cos_module.py  _cos_module.so*  cos_module_wrap.c  setup.py

现在我们能加载和执行cos_module，就好像我们之前做的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

In [1]: import cos_module

In [2]: cos_module?
Type:       module
String Form:<module 'cos_module' from 'cos_module.py'>
File:       /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/swig/cos_module.py
Docstring:  <no docstring>

In [3]: dir(cos_module)
Out[3]:
['__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 '_cos_module',
 '_newclass',
 '_object',
 '_swig_getattr',
 '_swig_property',
 '_swig_repr',
 '_swig_setattr',
 '_swig_setattr_nondynamic',
 'cos_func']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0

我们再次检验健壮性，看到得到了更好的错误信息(然而，严格地说Python中没有double类型)：

1
2
3
4
5
6
7

In [7]: cos_module.cos_func('foo')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-11bee483665d> in <module>()
----> 1 cos_module.cos_func('foo')

TypeError: in method 'cos_func', argument 1 of type 'double'

Numpy支持

numpy通过numpy.i文件提供了对SWIG的支持。这个接口文件定义了各种所谓的类型映射(typemaps)来转换Numpy数组和C数组。在下面的例子中我们将简略地看看这种类型映射在实际中如何起作用。

我们使用在ctypes例子中相同的cos_doubles函数：

1

void cos_doubles(double * in_array, double * out_array, int size);

1
2
3
4
5
6
7
8
9
10

#include <math.h>

/*  Compute the cosine of each element in in_array, storing the result in
 *  out_array. */
void cos_doubles(double * in_array, double * out_array, int size){
    int i;
    for(i=0;i<size;i++){
        out_array[i] = cos(in_array[i]);
    }
}

使用SWIG接口文件将它包裹为cos_doubles_func：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

/*  Example of wrapping a C function that takes a C double array as input using
 *  numpy typemaps for SWIG. */

%module cos_doubles
%{
    /* the resulting C file should be built as a python extension */
    #define SWIG_FILE_WITH_INIT
    /*  Includes the header in the wrapper code */
    #include "cos_doubles.h"
%}

/*  include the numpy typemaps */
%include "numpy.i"
/*  need this for correct module initialization */
%init %{
    import_array();
%}

/*  typemaps for the two arrays, the second will be modified in-place */
%apply (double* IN_ARRAY1, int DIM1) {(double * in_array, int size_in)}
%apply (double* INPLACE_ARRAY1, int DIM1) {(double * out_array, int size_out)}

/*  Wrapper for cos_doubles that massages the types */
%inline %{
    /*  takes as input two numpy arrays */
    void cos_doubles_func(double * in_array, int size_in, double * out_array, int size_out) {
        /*  calls the original funcion, providing only the size of the first */
        cos_doubles(in_array, out_array, size_in);
    }
%}

• 为了使用Numpy类型映射，需要numpy.i文件。
• 观察import_array()的调用，我们已经在Numpy-C-Api的例子中见到过。
• 因为类型映射仅仅支持参数ARRAY, SIZE我们需要包裹cos_doubles为cos_doubles_func，该函数接受两个数组包含各自大小作为输入。
• 相对于简单SWIG的例子，我们不需要包含cos_doubles.h头文件，因为我们通过cos_doubles_func暴露这个功能，我们没有其它东西想暴露给Python。

然后，如前述用distutils包裹它：

1
2
3
4
5
6

from distutils.core import setup, Extension
import numpy

setup(ext_modules=[Extension("_cos_doubles",
      sources=["cos_doubles.c", "cos_doubles.i"],
      include_dirs=[numpy.get_include()])])

显然，我们需要include_dirs指定位置。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

$ ls
cos_doubles.c  cos_doubles.h  cos_doubles.i  numpy.i  setup.py  test_cos_doubles.py
$ python setup.py build_ext -i
running build_ext
building '_cos_doubles' extension
swigging cos_doubles.i to cos_doubles_wrap.c
swig -python -o cos_doubles_wrap.c cos_doubles.i
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/esc/anaconda/include/python2.7 -c cos_doubles.c -o build/temp.linux-x86_64-2.7/cos_doubles.o
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/esc/anaconda/include/python2.7 -c cos_doubles_wrap.c -o build/temp.linux-x86_64-2.7/cos_doubles_wrap.o
In file included from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1722,
                 from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:15,
                 from cos_doubles_wrap.c:2706:
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION"
gcc -pthread -shared build/temp.linux-x86_64-2.7/cos_doubles.o build/temp.linux-x86_64-2.7/cos_doubles_wrap.o -L/home/esc/anaconda/lib -lpython2.7 -o /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/swig_numpy/_cos_doubles.so
$ ls
build/         cos_doubles.h  cos_doubles.py    cos_doubles_wrap.c  setup.py
cos_doubles.c  cos_doubles.i  _cos_doubles.so*  numpy.i             test_cos_doubles.py

接着，确信它起作用：

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pylab
import cos_doubles

x = np.arange(0, 2 * np.pi, 0.1)
y = np.empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab.plot(x, y)
pylab.show()

Cython

Cython不仅使用来写C扩展的Python样子的语言，而且是这个语言的一个高级编译器。Cython语言是Python的超集，包含额外的结构允许你调用C函数，将变量和类属性解释为C类型。在这个意义上可以叫它Python的一个类型。

除了这些几本的包裹原生代码的用例，Cython支持一个额外的用例，即交互优化。基本上是，从纯Python代码脚本出发逐步向代码瓶颈增加Cython类型来优化那些真正值得优化的代码。

在这个意义上它和SWIG非常相似，因为C代码可以自动生成，但某种意义上它也相当类似与ctypes，因为它包裹代码可以(几乎可以)用Python写成。

尽管其它自动生成代码方案会很难调试(例如SWIG)，Cython带有一个GNU调试器的扩展，能帮助调试Python，Cython和C代码。

注意：自动生成的C代码使用了Python-C-Api。

优势

• 类Python的语言来写C扩展
• 自动生成代码
• 支持增量优化
• 包含一个GNU调试器扩展
• 支持C++(自从0.13版本)

劣势

• 必须编译
• 需要额外的库(但仅仅在编译时，这个问题可以通过传递一个生成的C文件克服)

示例

我们cos_module的主要的Cython代码包含在文件cos_module.pyx中：

1
2
3
4
5
6
7

""" Example of wrapping cos function from math.h using Cython. """

cdef extern from "math.h":
    double cos(double arg)

def cos_func(arg):
    return cos(arg)

注意额外的关键字像cdef和extern。cos_func紧接着是纯Python。

我们再次使用标准distutils模块，但是这次我们需要一些来自Cython.Distutils额外的片段：

1
2
3
4
5
6
7

from distutils.core import setup, Extension
from Cython.Distutils import build_ext

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("cos_module", ["cos_module.pyx"])]
)

编译它：

1
2
3
4
5
6
7
8
9
10
11
12
13

$ cd advanced/interfacing_with_c/cython
$ ls
cos_module.pyx  setup.py
$ python setup.py build_ext --inplace
running build_ext
cythoning cos_module.pyx to cos_module.c
building 'cos_module' extension
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/include/python2.7 -c cos_module.c -o build/temp.linux-x86_64-2.7/cos_module.o
gcc -pthread -shared build/temp.linux-x86_64-2.7/cos_module.o -L/home/esc/anaconda/lib -lpython2.7 -o /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython/cos_module.so
$ ls
build/  cos_module.c  cos_module.pyx  cos_module.so*  setup.py

然后运行它：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

In [1]: import cos_module

In [2]: cos_module?
Type:       module
String Form:<module 'cos_module' from 'cos_module.so'>
File:       /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython/cos_module.so
Docstring:  <no docstring>

In [3]: dir(cos_module)
Out[3]:
['__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 '__test__',
 'cos_func']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0

接着，测试健壮性，可以看到我们获得了很棒的错误信息：

1
2
3
4
5
6
7
8
9

In [7]: cos_module.cos_func('foo')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-11bee483665d> in <module>()
----> 1 cos_module.cos_func('foo')

/home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython/cos_module.so in cos_module.cos_func (cos_module.c:506)()

TypeError: a float is required

另外，值得注意的是Cython带有完整的C数学库声明，将之上代码简化为：

1
2
3
4
5
6

""" Simpler example of wrapping cos function from math.h using Cython. """

from libc.math cimport cos

def cos_func(arg):
    return cos(arg)

在这个例子中cimport声明被用来importcos函数。

Numpy支持

Cython通过numpy.pyx文件支持Numpy，这允许你将Numpy数组类型添加到Cython代码。例如将i指定为int类型，将变量a指定为numpy.ndarray并给定dtype。某些优化像边界检查也支持。参看Cython文档的相关章节。万一你想将Numpy数组作为C数组传递给你的Cython包裹的C代码，Cython维基中有一个章节。

在以下例子中，我们将展示如何如何使用Cython包裹熟悉的cos_doubles函数。

1

void cos_doubles(double * in_array, double * out_array, int size);

1
2
3
4
5
6
7
8
9
10

#include <math.h>

/*  Compute the cosine of each element in in_array, storing the result in
 *  out_array. */
void cos_doubles(double * in_array, double * out_array, int size){
    int i;
    for(i=0;i<size;i++){
        out_array[i] = cos(in_array[i]);
    }
}

该函数使用以下Cython代码被包裹为cos_doubles_func：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

""" Example of wrapping a C function that takes C double arrays as input using
    the Numpy declarations from Cython """

# import both numpy and the Cython declarations for numpy
import numpy as np
cimport numpy as np

# if you want to use the Numpy-C-API from Cython
# (not strictly necessary for this example)
np.import_array()

# cdefine the signature of our c function
cdef extern from "cos_doubles.h":
    void cos_doubles (double * in_array, double * out_array, int size)

# create the wrapper code, with numpy type annotations
def cos_doubles_func(np.ndarray[double, ndim=1, mode="c"] in_array not None,
                     np.ndarray[double, ndim=1, mode="c"] out_array not None):
    cos_doubles(<double*> np.PyArray_DATA(in_array),
                <double*> np.PyArray_DATA(out_array),
                in_array.shape[0])

可以使用distutils编译：

1
2
3
4
5
6
7
8
9
10

from distutils.core import setup, Extension
import numpy
from Cython.Distutils import build_ext

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("cos_doubles",
                 sources=["_cos_doubles.pyx", "cos_doubles.c"],
                 include_dirs=[numpy.get_include()])],
)

• 如上编译Numpy的例子，我们需要include_dirs选项。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

$ ls
cos_doubles.c  cos_doubles.h  _cos_doubles.pyx  setup.py  test_cos_doubles.py
$ python setup.py build_ext -i
running build_ext
cythoning _cos_doubles.pyx to _cos_doubles.c
building 'cos_doubles' extension
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/esc/anaconda/include/python2.7 -c _cos_doubles.c -o build/temp.linux-x86_64-2.7/_cos_doubles.o
In file included from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1722,
                 from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17,
                 from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:15,
                 from _cos_doubles.c:253:
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION"
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/__ufunc_api.h:236: warning: ‘_import_umath’ defined but not used
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include -I/home/esc/anaconda/include/python2.7 -c cos_doubles.c -o build/temp.linux-x86_64-2.7/cos_doubles.o
gcc -pthread -shared build/temp.linux-x86_64-2.7/_cos_doubles.o build/temp.linux-x86_64-2.7/cos_doubles.o -L/home/esc/anaconda/lib -lpython2.7 -o /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython_numpy/cos_doubles.so
$ ls
build/  _cos_doubles.c  cos_doubles.c  cos_doubles.h  _cos_doubles.pyx  cos_doubles.so*  setup.py  test_cos_doubles.py

如前述确保它能起作用：

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pylab
import cos_doubles

x = np.arange(0, 2 * np.pi, 0.1)
y = np.empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab.plot(x, y)
pylab.show()

总结

这个章节中四种不同和本地代码接口技术被呈现在您面前。这个表格简要的总结了这些技术的某些方面。

x |Part of CPython |Compiled |Autogenerated |Numpy Support
|:——|:—-|:—–|:—–|:—-|
Python-C-Api |True |True |False |True |
Ctypes |True |False |False |True |
Swig |False |True |True |True |
Cython |False |True |True |True |

相比所有技术中，Cython是最现代最高级的了。特别是，通过向Python代码中添加类型增量优化代码的能力是独一无二的。

更多阅读和参考

• Gaël Varoquaux的关于避免复制数据的博文对于如何聪明的处理内存管理提供了一些深刻的见解。如果你曾经碰到过大数据集的问题，这是一个可以激发灵感的参考。

练习

因这是一个全新的章节，这些练习更可视为是接下来看什么的指针。所以选择您最感兴趣的那个。如果你对此有更多好的想法，请联系我们！

1. 下载每个练习的源码并在你的机器上编译运行。
2. 对每个例子做小的修改确信它能起作用。(像将cos改成sin)

3. 大多数例子，特别是涉及Numpy的例子可能仍然对错误输入很脆弱、并且返回模糊的消息。寻找使这些例子出问题的方法，指出问题是什么并且设计潜在的解决方案。这有一些提示：

   1. 数值溢出
   2. 输入输出数组长度不同
   3. 多维数组
   4. 空数组
   5. 非double型数组

4. 使用IPython中的magic%timeit来测量不同方案的执行时间。

Python-C-Api

1. 更改Numpy的例子让函数接受两个输入参数，第二个参数是预分配的输出数组，让它像其它Numpy例子。
2. 更改例子让函数仅仅接受一个输入数组并且原地修改。
3. 尝试使用Numpy迭代协议修正例子，如果你设法获取了一个工作的解决方案，请在github上发布一个拉取请求(pull-request)。
4. 你也许注意到了，Numpy-C-API例子是唯一不包裹cos_double但是直接将cos应用到Numpy数组的元素的Numpy例子。这相对于其它技术有什么优势？
5. 你能仅仅使用NJumpy-C-API包裹cos_doubles吗？你可能需要确保数组是正确的类型，并且在内存中一维连续。

Ctypes

1. 更改像cos_double_func的Numpy例子为你处理预分配，使之更像Numpy-C-Api例子。

SWIG

1. 查看SWIG自动生成的代码，你能理解多少？
2. 更改Numpy的例子像cos_double_func处理预分配，让它更像Numpy-C-API的例子。
3. 更改C函数cos_doubles让它返回一个分配的数组。你能用SWIG类型映射包裹它？如果不能，为何不行？有没有特定条件的变通方案。(提示：你知道输出数组的大小，所以可能从返回的double *构建一个Numpy数组。)

Cython

1. 查看Cython自动生成的代码。仔细看看Cython插入的一些注释。你看到了什么？
2. 查看Cython文档的章节Working with Numpy去学习使用Numpy如何增量优化纯python代码。
3. 更改Numpy例子比如cos_doubles_func处理预分配，使之更像Numpy-C-Api例子。

FootNotes

[^1]:我不懂，还没看Advanced Numpy部分

不要相信我说了什么，我所说的 未必 是事实，你也 未必 能理解成事实。

总之，经历过回头看就自知了。

莫听穿林打叶声，何妨吟啸且徐行。竹杖芒鞋轻胜马，谁怕，一蓑烟雨任平生。

料峭春风吹酒醒，微冷，山头斜照却相迎。回首向来萧瑟处，归去，也无风雨也无晴。

苏轼，定风波

目录

• toc
  {: toc}

初试

我不记得为什么会做出这些选择，如果非要找个理由。正如Steve Jobs2005年在斯坦福大学演讲时所说的，I follow my heart.

缘起

因为那本叫《Debian标准教程》的书作者曾是byr linux版版主？

看过一本讲TeX的书的作者提到了byr？

然后我来到了byr linux版……

然后接到母亲的电话说可以考北邮，北邮的光学好……我说好啊……

然后准备选通信原理作为专业课，因为我想跨考通信。

买了本书，看了两个月不知道什么玩意……

论坛上有人说，喜欢编程考信安吧，考804信号与系统简单。复试无所谓，大家都不会。

然后不知道为什么我就做了这辈子可能都无法理解的选择……我好像喜欢的不是编程……

复习

八月初说准备考研吧。买本新东方乱序版单词，结果背了一半不到……然后，没有然后了，忘干净了。

开始准备复习高数，李永乐全书还是从学长那要的去年的。这么断续看了一个月，被虐得不轻。好像题目就没有会做且做对的。最后我只看例题了……

专业课一两个月看完了，于是扔掉了……

买本政治吧，大纲解析真贵。我觉得应该最后一个月再看吧，没想到最后根本就没看……

风中劲草买的较早吧，不去自习每天就看看劲草，在寝室读着玩。进度慢爆，拖拖拉拉，同学都复习两遍三遍了我第一遍都没。

打酱油一样的自习，早上八九点起，十一点回，下午三点多去，五点回，晚上不去……一路游击，占个毛线座。

好不容易刻苦了两天，然后生病了。紧接着一个半月书都没碰，天天自己逛公园和看电影去了。

最后只二十来天，想想数学真题看看吧，全书我就看了看例题。历年真题，基本上没答对几道，基本上不是出错就是根本不会……

最后两周，忽然想起貌似英语还没看……把早就买了的真题黄皮书以一天50-180页的速度看完……然后花大概七天背了背附送的单词书，不过貌似没啥用。

考试

运气很好分到个有空调教室，还不靠门。并不是所有人运气都这么好，没空调分阴面的同学在无锡湿冷湿冷的冬天冻爆了。

考前晚上还忙着为恶搞百度百科P图。一副三心二意的模样，考前一周竟然还忙着折腾scipy……

考政治那天晚上漫天大雪，想到要考完试了高兴异常。我又到公园转了一圈才去自习，没怎么看就回寝室，读了十分之一十八大报告就跪了，胡总真不容易。

晚上睡觉时楼上吵得要死，很烦的想带上耳塞等会去掉，结果一戴上睡着了……

第二天下雪了，政治做得没什么特殊感觉，有点冷。

英语，基本都不会，一路瞎蒙，作文划的一塌糊涂。

第二天兴奋的不知两点三点才睡着，想着要结束了就高兴。

上午数学做得史无前例的顺，做完了看还有半个多小时，又翻翻看看，忽然发现有误解的地方，胡乱一改然后就收卷了。然后听旁边同学说数学难……

专业课前半小时后半小时净听监考老师blabla个没完没了，然后简单的会做难得一点看不懂。感觉有的题还是错的。

亢奋的情绪持续了第二天一天，晚上又是两三点才睡着。高兴过头了。

复试完第一件事竟然是装gentoo……太TM屌丝了……

复试

十一号报到十号坐高铁滚到北京，可是奢侈了一把。承蒙同学收留，住在了北科。复试这几天爽爆了，日行二三十公里……北科距北邮40分钟步程(我的走路速度就和慢跑差不多了……)一天至少也得往来一趟。在北邮瞎转一圈，好在北邮小到能从南门看到北门西门看到东门，没觉得太累。第二天又去中科院那边找同学，一路问一路走饶了一大圈又爽了一把，晚上又和同学逛到九点多……但谢谢他(她)们肯花时间招待我，陪我瞎转，很温暖的北京之旅啊。明天打算去北大和清华混吃求陪去，然后哈哈，亦庄。

体检

早上sb了动作慢了一步，同寝的国防生们早冲出去跑操去了。结果面对北科大宿舍新装的伟大门禁我就sb了……等到一会有收垃圾的开了个侧门我终于溜出宿舍了……

结果就是，本打算6点多去校医院排队体检，结果七点才离开北科，又花了四十分钟到北邮校医院。倒是远远看到一长对一路排到教三。然后竟然好多排队不拿表……

拿表，填好，听说可以先体检后缴费就四处找人少的地方。结果后来就sb了，差点没测完……

建议是早去先胸透视力然后抽血，男外别去太晚，其它好像都特快。不然就没人给您体检了……说到男外科我要扯扯，人说北邮校医院有一大绝技——切包皮。我感觉男外那男医生净督促大家切包皮了。

笔试

唯一印象深刻的是我爸发短信说去年要带2b铅笔，让我带上。我根本就不知道，差点悲剧掉。

据说信安的笔试是浮云，我也无心准备。看了看计网，两天内忘光。看看c语言，忘得更快，好歹毕设还和C有些关系。准备看概率论，最后一个字没看……

笔试一个小时，五选三，90道选择题，都是涂答题卡。只要熟悉都不难。

官方说明复试挺重要，由导师自主划线，决定复试人员。但没看到有任何人因为笔试挂掉。

反正没啥好担心的，把复习的时间省下来做其它有意义的事吧。

唯一突然的是，笔试结束后留下来填张表，关于考研成绩，本科专业，是否服从分配，认同团队文化，最早什么时候能来，是否同意调剂到专硕等问题。钮老师还提醒了下复试要带上缴费后打印的报名表，成绩单和简历和刚发的复试表。

面试

面试下午一点半开始，我成绩单简历报名表什么都没打，还放在北科。于是花了四十分钟又从北邮回到北科。半路接到辅导员电话，说查不到入党时间，要问之前负责的C老师，老师听说我在帝都复试，很和善的祝愿我了，瞬间感觉很温暖……回到北科进不去寝室，躺在操场的长椅上看国防生训练，然后就睡着了……

醒来想起入党时间还不知道，打电话问C老师，结果基本是被骂了一遍，呼呼拉拉又打了几通电话，最后老师说，好吧，下午四点半之后打给我吧。结果下午复试时收到老师短信告诉我是几号，真假不知，不过觉得不会坑我，我也不在乎，总之，谢谢C老师。

摸摸兜里就六七块钱了，在北科随便吃了顿饭。然后终于能回寝室了，打印好东西直接又花了四十分钟回到北邮，好在gxy之前带我在北邮逛了一圈，没问几个人就找到新科研楼12层了。

到上面发现我是第二个……qq群里问心是第一个，他进去我就在外面跟一起等待的人扯淡开玩笑。然后我就进去了。

不是什么群面，单独面试感觉挺好。对面一溜五六个老师，我对面坐下。

挺囧的面试……开始让一分钟英文自我介绍，题材随意。我讲了几句看没人听。于是说老师不会的单词用中文行么。老师说算了，你别说了……

然后就开始问我的特长，对什么感兴趣，哪方面研究深入。看我简历说熟悉linux，问我什么是驱动？我完全没明白老师想问什么，很坦然说不知道。出去才想到好像是问我驱动做了什么……也许老师感觉他问了个linux很基本的问题，我都不会，一定在瞎忽悠吧。

发的文章怎么回事，做了多少事，是不是第一作者。项目中你做了多少做了什么……

兴趣点太多，简历写得跳跃度太大。被直接看出来我是很冲动型的人了。不过说我什么坚持不太久我就不高兴了，单linux用这么久了都（好吧，快两年而已）。

问问我成绩怎样，班级排名怎样。班里考研的情况怎样。

问了问我用没用过gcc，用linux能不用gcc？又问用没用过gdb，怎么用，我随口说用过，是的，怎么用gdb调试python我知道，但从来没用gdb调试过c,我说在python中用过(翻译scipy lecture notes时学的)……但老师肯定不知道我的意思……总之，交流能力当时sb了。

然后扯到毕设，好像我扯过去的……扯扯貌似大家没怎么明白，就没怎么问。

问了问为什么选信安，好吧，我就是选信安而已，没啥理由。面对这神回答，老师很无语，我也很无语……

问了我觉得自己有啥优缺点。我说那动手能力和学习能力算有点吧，缺点，有点胆小吧。老师说我面了这么多学生没感觉你胆小，我说我心理有准备。又问问有没其它缺点，我想想说脾气不太好吧，老师说和同学关系不好？我说没有啊，只是有时候感觉太关注自己忽视别人了……

又问如果碰到分到不喜欢的工作会不会很别扭，你这看上去太冲动型了。我说老师您说得对，但我也没啥不感兴趣的，两年半总坚持的了。又说毕竟这个年龄想多看看多学学多折腾折腾，老师点点头，也不知道啥意思。

问我有啥不顺心的事吧，我实在想不到学校有啥，结果就问家里，结果我说了什么来着最后老师竟然来句对不起，我妹的我sb的回了句没事……

问了我到底对啥比较感兴趣，我说数学和计算机吧，渗透性比较强，我也喜欢各种各样交叉的东西，比如说信安。这时候说也有可能让你写文档，我说，好啊，老师您看简历，我很喜欢写东西啊，像我用的archlinux的社区维基上我就写了和翻译了很多东西啊。

被诘难是否没做出什么东西，我说老师您看简历上，做了很多小玩意，虽然不是啥特别好的东西。好歹是我一点一点查资料尝试错误，一点一点实现的。

然后走得时候老师们都很和善，然后我很sb的溜走了，出门被一群人围上了，怎么你这么久，好吧，我感觉就一会的事。

然后……想起来包忘里头了，手里还拿着不知道谁的笔……得，一边跟别人闲扯复试很简单，淡定加油或者开玩笑，一边等下一位出来。下一个出来时很猥琐地进去把包拿走笔放下了……我不知道老师怎么想……这孩子太sb了……

好了，这就是我的研究生复试面试了。sb的人生不需要解释,猥琐的举止无需掩饰……

政面

一早赶去没什么人，老师看看表填好没，随便问问个人信息性格政治倾向一类的。很快。

总结

• 首先，做什么都要有信心吧。不要总是消极自我暗示，要积极自我欺骗。
• 多活动活动，劳逸结合，效率比时间重要。时间是像手中的沙子，抓得越紧，从指缝间流失的越快。
• 考研没什么重要的，身康体健，亲友健在，才是最重要的。

• 别跟大师比较，人比人，比死人。仅仅如一首禅诗：

  To follow the path:（沿着这样一条道路：） look to the master,（寻找大师，） follow the master,（跟随大师，） walk with the master,（与大师通行，） see through the master,（洞察大师，） become the master.（成为大师。）

• 尽人事，听天命。不要胆怯，如果你真的在乎。

• 拿得起，放的下。敢赌敢输，天下之大，何处无容身之地。
• 福兮祸之所倚,祸兮福之所伏。

致谢

• 我的亲人：我不在乎上不上研，只在乎你们。
• 北邮人考研专版：很棒的地方。
• 同学及朋友：谢谢鼓励，谢谢复试时的款待、收留和帮助。
• 老师：从未曾为难过我

总之，很高兴，某件事又结束了。

二维相位去包裹：路径积分法

本文介绍几种基于路径积分的二维相位展开方法。这只是个预告，表示这几天没有闲着。

##为什么二维相位展开会发生问题？

Itoh’s算法的二维扩展

Itoh’s一维相位展开算法很容易扩展到二维。

我们可以先对第一列进行解包裹，然后在此基础上对每行进行解包裹，整个相位图像就都能被解包裹了。

我把扩展的itoh’s算法函数写在util中并导入

1
2

from util import *
from mpl_toolkits.mplot3d import axes3d

获取实验相位函数。我随便选的，应该没问题吧……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

X, Y, Z = axes3d.get_test_data(0.05)
Z_origin = Z.copy()
plt.figsize(10,10)

plt.subplot(121)
plt.title("origin phase")
plt.imshow(Z, cmap='gray')
ax = plt.gca()
ax.set_axis_off()
colorbar(shrink=.40)

plt.subplot(122)
plt.title("wraped phase")
ax = plt.gca()
ax.set_axis_off()
Z = wrap(Z)
plt.imshow(Z, cmap='gray')
colorbar(shrink=.40)

紧接着用itoh的方法解包裹

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

Z = raster_unwrap2(Z)

plt.figsize(10,10)

plt.subplot(121)
ax = plt.gca()
ax.set_axis_off()
plt.title("unwraped phase")
plt.imshow(Z, cmap='gray')
colorbar(shrink=.40)

plt.subplot(122)
ax = plt.gca()
ax.set_axis_off()
plt.title("erro map")
plt.imshow(Z - Z_origin, cmap='gray')
colorbar(shrink=.40, ticks=None)
np.all(np.abs((Z - Z_origin) < 1e-8))

正如我们所见，这算法效果很好。但是如果数据点出了点问题的话……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

Z = wrap(Z)

x = np.random.randint(0,120,100)
y = np.random.randint(0,120,100)
Z[x,y] = random.random() * 3

plt.subplot(131)
ax = plt.gca()
ax.set_axis_off()
plt.title("noisy wraped phase")
plt.imshow(Z, cmap='gray')
colorbar(shrink=.25)

Z = raster_unwrap2(Z)

plt.subplot(132)
plt.title("noisy unwraped phase")
ax = plt.gca()
ax.set_axis_off()
plt.imshow(Z, cmap='gray')
colorbar(shrink=.25)

np.all((Z == Z_origin))

plt.subplot(133)
plt.imshow(Z-Z_origin, cmap='gray')
ax = plt.gca()
ax.set_axis_off()
plt.title("error map")
colorbar(shrink=.25)

我加了一些随机噪声，结果就悲剧了。查看误差图给我们了一些启示。在某些出问题的点同一行又边全都出了问题。这种解包裹误差不断扩大并形成一条线的现象，似乎 叫做 拉线 现象(我瞎说的，只是有文献提到过拉线现象，也可能就是这个吧……)。

如果我们对Itoh算法换一种方式扩展。先对第一行解包裹，再在此基础上对各个列解包裹。如果没有有问题的点。结果和第一种扩展方式一样。这给我们个启示，我们可以在任意方向进行相位去包裹，但只要碰到出问题的点都会在之后的积分路径中出现问题。

如果要解包裹出正确的相位，必须绕过这些奇点(有问题的点)。

事实上，相位解包裹领域借用了高等数学和复变函数分析里的理论和观点，并最终形成了相位解包裹问题的留数定理。

相位去包裹留数定理

要想解包裹与路径无关，必须有：

$$\oint \nabla \varphi(r) \cdot dr = 2\pi \times \text{(sum of enclosed residue charges)}$$

在数字图像中，最小的环路是2X2大小的四个像素，如下图所示。通过计算$\sum_{i=1}^4\Delta_i$是否为零，为零则无留数，否则此区域内有正或负留数。实际上，这些有问题的点留数通常不会为0。

为了在任意方向上能进行正确的解包裹，必须想办法使环路内的留数和为0,让积分路径(也是解包裹的路径)绕过这些留数不为0的奇点。

之后所提到的，我所看到的一些路径积分方法。都是通过尽量在可靠的地方进行相位展开(即积分),想办法绕开可能造成问题的点，或平衡这些造成问题的点。

之后做好准备看看那些大牛们发展出的各种路径积分算法吧。我已经佩服的五体投地了。这复杂程度已经让个人Keep it simple， stupid的人生哲学彻底崩坏。各路人马的算法涉及大量数学理论、信号处理理论(mask、noise filter、quality map)、图形形态学操作(开闭操作、连通性判断)、图论(分支切割算法)、最小生成树……完全看跪了……如果你准备动手实践，为了提高时间和空间效率，还有涉及可能是算法设计的东西吧[^1]。

几种路径积分算法

Goldstein 分支切割(branch cut)算法

该算法旨在用最短的路径平衡正负留数，设置解包裹时的障碍，使解包裹路径能按照正确的方向进行。总体上说，该算法相当高效，很容易查看效果。但跟质量图没关系，很容易错误的放置障碍，导致错误的积分路径。

质量引导(quality-guided)算法

先积分质量高的地方，后积分质量低的地方，这基本能保证解包裹路径的正确。这个算法和留数没关系，所以不能保证不出错。但如果有个非常好的质量图，结果相当棒。另外，效率也不错。

掩码分支(mask)算法

前两种算法的杂交，通过某种方法沿着质量最低的地方扩展mask，最后解包裹时绕过这些mask。同上方法，没有好的质量图就别去试。

flynn最小断点算法

很直观的算法。就是把不连续的部分通过增加整数倍个$2\pi$变得连续，虽然说起来没这么简单。该算法用了个生成树，效果很好，还可以用质量图获得更好的效果……

Off-Topic

相位去包裹笔记放在这里。

Footnotes

[^1]:在 Data Analysis with Open Source Tools 一书中推荐过 一本书叫如何设计算法 ，应该值得一看。

一维相位去包裹：原理与仿真

Off-Topic

未觉池塘春草色，阶前梧叶已秋声。

——朱熹

时光飞逝，转眼就要毕业了。前天导师打电话摧我，这几个月就毕业了，你那毕业设计得赶紧做了。想想真是，四年倏忽过去，浑然不觉，已近毕业。

也许这是最后自在平静的日子了，时间却像沙子，越用力抓紧，从指间溜走地越快。

两年过得太快，I follow my passion，做了太多没用的事，有得有失。只是希望：

You can’t connect the dots looking forward you can only connect them looking backwards. So you have to trust that the dots will somehow connect in your future. You have to trust in something: your gut, destiny, life, karma, whatever. Because believing that the dots will connect down the road will give you the confidence to follow your heart, even when it leads you off the well worn path.

——Steve Jobs, Stanford Commencement Adress, 2005

为何相位重要

我大学的专业是光学，时至今日，依然清晰记得两年前在阴暗的实验室做全息摄像的实验。我们几个同学花了好久把光路摆好，让激光打在物体上反光到干版上曝光。然后拿到什么试剂中定形。最后在拿出来用激光将物体的影响重新清晰的放出三维的像时，心里相当兴奋。

为什么全息照相能产生立体的像呢？普通的照相技术都是仅仅记录光波的强度，而不记录相位，因此失去了很多相位中的信息。但全息照相通过相干光(激光)之间的干涉在干版上同时记录下了强度和相位，再用相干光照射干版重放，好像光是从真的物体发出的一样。

如果您也做过这种实验，您应该知道为什么相位重要了。

这还有另一个例子，关于爱因斯坦和蒙娜丽莎。

我从《Two Dimensional Phase Unwrapping Theory Algorithms and Software》中看到了这个例子，然后自己动手用python试了试。

以下代码用来交换两个图片的相位

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

from scipy.fftpack import *

# <codecell>

# read file
im_1 = plt.imread('einstein.jpg')
im_2 = plt.imread('monalisa.jpg')

# fft and reverse two images' phase
m_1, p_1 = np.abs(fft2(im_1)), np.angle(fft2(im_1))
m_2, p_2 = np.abs(fft2(im_2)), np.angle(fft2(im_2))

# <codecell>

im_swapphase_1 = np.real(ifft2(m_1 * np.cos(p_2) + m_1 * np.sin(p_2) * 1j))
im_swapphase_2 = np.real(ifft2(m_2 * np.cos(p_1) + m_2 * np.sin(p_1) * 1j))

# <codecell>

plt.figsize(10,10)
plt.gray()
plt.subplot(2,2,1)
plt.imshow(im_1, origin='lower')
plt.subplot(2,2,2)
plt.imshow(im_2, origin='lower')
plt.subplot(2,2,3)
plt.imshow(im_swapphase_1, origin='lower')
plt.subplot(2,2,4)
plt.imshow(im_swapphase_2, origin='lower')

上面两幅图是交换相位前的图，下面两幅是之后的。显然，一团糟，相位中是有信息的。

为何要相位去包裹

简单地说，就是说，任何仪器，比如说量角器，顶多只能测得($-\pi, \pi$]之间的量,但真正的相位角度不该这样，而是分布在实数空间内，应该是测得的($-\pi, \pi$]之间的值的2$\pi$整数倍。

真正的相位值被“包裹”起来了，但为什么要解包裹呢？

不解包裹无法对相位进行计算。

Itoh的路径积分法

我们先定义一个获取包裹相位的算子$\mathscr{W}$，该算子将相位包裹，获取位于$(-\pi,\pi]$之间的包裹相位。

$$\mathscr{W}\varphi = \arctan[\cos(Real \varphi / Img \varphi)]$$

还可以这样写

$$\mathscr{W}{\varphi(n)} = \psi(n) = \varphi(n) + 2\pi k(n)$$

其中$k(n)$是使包裹相位位于$-\pi,\pi$之间的值。

显然包裹相位$\psi(n)$有：

$$\pi \geq \psi(n) \gt -\pi$$

定义差分算子$\Delta$：

$$\Delta {\varphi(n)} = \varphi(n+1) - \varphi(n)$$

$$\Delta {k(n)} = k(n+1) - k(n)$$

计算被包裹相位的差分：

$$\Delta {\psi(n)} = \Delta {\varphi(n)} + 2\pi\Delta {k_1(n)}$$

我们再用$\mathscr{W}$作用于该差分得：

$$\mathscr{W}{\Delta {\psi(n)}} = \Delta {\varphi(n)} + 2\pi[\Delta {k_1(n)} + k_2(n)]$$

显然上式结果应该位于$(-\pi,\pi]$，假如此时还有$\Delta {\varphi(n)}$也位于$(-\pi,\pi]$,则上式右边第二项$2\pi[\Delta {k_1(n) } + k_2(n)]$应该为零，则有：

$$\Delta {\varphi(n)} = \mathscr{W}{\Delta {\psi(n)}} $$

显然，由该差分式可得：

$$\varphi(m) = \varphi(0) + \sum_{n=0}^{m-1} \mathscr{W}{\Delta {\mathscr{W}{\varphi(n)}}} $$

上式说明，真实相位可以通过对包裹相位的差分的包裹进行积分求得。

于是itoh的一维相位去包裹算法综述如下：

对信号相位数组$\psi(i),0 \leq i \leq N-1$

• 计算相位差分$D(i) = \psi(i+1)-\psi(i), i=0,\ldots,N-2$
• 计算包裹的相位差分$\Delta(i) = \Delta{D(i)}, i=0,\ldots,N-2$
• 初始化初值$\varphi(0)= \psi(0)$
• 累加解包裹$\varphi(i) = \psi(i) + \Delta(i)$

Itoh的方法很简单实用，但受到两个重要因素的影响：相位失真和噪声。下面仿真两种情况的影响。

仿真

对正弦波相位函数(间谐波，一切波的基础)

$$\varphi(t) = 10\sin(10t), 0 \leq t \leq 1$$

通过计算可以得知，使之不产生相位失真，区间内至少有32个采样点:

对相位变化有

$$\Delta \varphi = \dot{\varphi}\Delta t$$

其中 $\dot{\varphi} = d\varphi / dt = 100\cos10t$ ，可看出相位在 $n\pi/10$ 取得极值，加上条件：

$$\left\vert\Delta \varphi\right\vert \lt \pi$$

则

$$\left\vert \frac{100}{N} \cos 10t \right\vert \lt \pi$$

得

$$N \gt 31.83$$

则对这个函数如果采样率低于32就会产生相位失真。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

import numpy as np
import matplotlib.pyplot as plt

# <codecell>

# phase function
def pf(x):
    return 10 * np.sin(10 * x)
# wrap function
def wrap(x):
    return np.arctan2(np.sin(x), np.cos(x))
# unwrap function
def wrap_diff(x):
     return wrap(np.diff(x))

def unwrap(x):
    y = x
    y[0] = x[0]
    for i in range(len(x) - 1):
        i += 1
        y[i] = y[i - 1] + wrap_diff(x)[i - 1]
    return np.array(y)

def noise(x, snr):
    return x + np.random.normal(loc=0.0, scale = np.sqrt(np.max(x) / 2 ** snr), size=len(x))

def show_unwrap(x,y,y_w,t,p):
    plt.ylim((-12,12))
    plt.ylabel("Phase in Radians")
    p_w, = plt.plot(x, y_w,'o:')
    p_o, = plt.plot(t,p,':')
    p_u, = plt.plot(x,y,'s')
    plt.legend([p_w, p_o, p_u], ["sampled wrapped phase", "original phase function", "unwrapped phase"])
    

# <headingcell level=1>

# 当取样点为num时

# <codecell>

# Origin
t = np.arange(0,1,0.01)
p = pf(t)

# plot
## num 50
plt.subplot(311)
plt.title("num = 50")
# sampled data
num = 50
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
## num = 32
plt.subplot(312)
plt.title("num = 32")
# sampled data
num = 32
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
## num = 31
plt.subplot(313)
plt.title("num = 31")
# sampled data
num = 31
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
plt.xlabel("Relative Time")


# <headingcell level=1>

# 噪声

# <codecell>

# noise influence
x = np.linspace(0,1,200)
y = pf(x)
y_10 = unwrap(wrap(noise(y,10))) + 20
y_5 = unwrap(noise(wrap(y), 5)) + 40
y_2 = unwrap(noise(wrap(y), 2)) + 60
y_1 = unwrap(noise(wrap(y), 1)) + 80
# plot
plt.xlabel("Relative Time")
plt.ylabel("Phase in Radians")
p_o, = plt.plot(t,p)
p_1, = plt.plot(x,y_1,'--')
p_2, = plt.plot(x,y_2,'-.')
p_5, = plt.plot(x,y_5,':')
p_10, = plt.plot(x,y_10,':')
plt.legend([p_o, p_1, p_2, p_5, p_10], ["Origin", "SNR=1", "SNR=2", "SNR=5", "SNR=10"])

# <codecell>

我们可以看到：

• 在采样率低于某一关键值32时基本没法解出正确相位。

  

• 信噪比越低，解包裹效果越差。
  

Summary

如果是二维相位去包裹问题，还有个奇点的问题。噪音、相位失真、奇点，似乎是所有相位解缠算法必须面对的三大问题。

相关Ipython Notebook文件和Einstein和Monalisa的图像在此下载

高级Python结构

原谅渣翻译，可能仅仅是给我自己看的。本来多年(也就几个月吧)之前将此文投递到OSChina翻译频道结果被以代码太多为由被拒，于是译者只好用自己的渣英语渣水平翻译给自己看了……，期待各路大婶们指正……

翻译自http://scipy-lectures.github.com/advanced/advanced_python/index.html

作者: Zbigniew Jędrzejewski-Szmek

这章有关Python中被认为高级的特性——就是说并不是每个语言都有的，也是说它们可能在更复杂的程序或库中更有用，但不是说特别特殊或特别复杂。

强调这点很重要：这一章仅仅关于语言自身——关于辅之以Python的标准库功能的特殊语法所支持的特性，不包括那些智能的外部模块实现。

在开发Python程序语言的过程中，它的语法，独一无二。因为它非常透明。建议的更改通过不同的角度评估并在公开邮件列表讨论，最终决定考虑到假设用例的重要性、添加更多特性的负担，其余语法的一致性、是否建议的变种易于读写和理解之间的平衡。这个过程由Python Enhancement Proposals(PEPs)的形式规范。最终这一章节中描述的特性在证明它们确实解决实际问题并且使用起来尽可能简单后被添加。

目录

• toc
  {: toc}

迭代器(Iterators), 生成表达式(generator expressions)和生成器(generators)

迭代器

简单

重复工作是浪费，将不同“土生土长”的方法替换为标准特性换来的是更加易于阅读和操作。

Guido van Rossum — Adding Optional Static Typing to Python

迭代器是依附于迭代协议的对象——基本意味它有一个next方法(method)，当调用时，返回序列中的下一个项目。当无项目可返回时，引发(raise)StopIteration异常。

迭代对象允许一次循环。它保留单次迭代的状态(位置)，或从另一个角度讲，每次循环序列都需要一个迭代对象。这意味我们可以同时迭代同一个序列不只一次。将迭代逻辑和序列分离使我们有更多的迭代方式。

调用一个容器(container)的__iter__方法创建迭代对象是掌握迭代器最直接的方式。iter函数为我们节约一些按键。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

>>> nums = [1,2,3]      # note that ... varies: these are different objects
>>> iter(nums)                           
<listiterator object at ...>
>>> nums.__iter__()                      
<listiterator object at ...>
>>> nums.__reversed__()                  
<listreverseiterator object at ...>

>>> it = iter(nums)
>>> next(it)            # next(obj) simply calls obj.next()
1
>>> it.next()
2
>>> next(it)
3
>>> next(it)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration

当在循环中使用时，StopIteration被接受并停止循环。但通过显式引发(invocation)，我们看到一旦迭代器元素被耗尽，存取它将引发异常。

使用for...in循环也使用__iter__方法。这允许我们透明地开始对一个序列迭代。但是如果我们已经有一个迭代器，我们想在for循环中能同样地使用它们。为了实现这点，迭代器除了next还有一个方法__iter__来返回迭代器自身(self)。

Python中对迭代器的支持无处不在：标准库中的所有序列和无序容器都支持。这个概念也被拓展到其它东西：例如file对象支持行的迭代。

1
2
3

>>> f = open('/etc/fstab')
>>> f is f.__iter__()
True

file自身就是迭代器，它的__iter__方法并不创建一个单独的对象：仅仅单线程的顺序读取被允许。

生成表达式

第二种创建迭代对象的方式是通过 生成表达式(generator expression) ，列表推导(list comprehension)的基础。为了增加清晰度，生成表达式总是封装在括号或表达式中。如果使用圆括号，则创建了一个生成迭代器(generator iterator)。如果是方括号，这一过程被‘短路’我们获得一个列表list。

1
2
3
4
5
6

>>> (i for i in nums)                    
<generator object <genexpr> at 0x...>
>>> [i for i in nums]
[1, 2, 3]
>>> list(i for i in nums)
[1, 2, 3]

在Python 2.7和 3.x中列表表达式语法被扩展到 字典和集合表达式。一个集合set当生成表达式是被大括号封装时被创建。一个字典dict在表达式包含key:value形式的键值对时被创建：

1
2
3
4

>>> {i for i in range(3)}   
set([0, 1, 2])
>>> {i:i**2 for i in range(3)}   
{0: 0, 1: 1, 2: 4}

如果您不幸身陷古老的Python版本中，这个语法有点糟：

1
2
3
4

>>> set(i for i in 'abc')
set(['a', 'c', 'b'])
>>> dict((i, ord(i)) for i in 'abc')
{'a': 97, 'c': 99, 'b': 98}

生成表达式相当简单，不用多说。只有一个陷阱值得提及：在版本小于3的Python中索引变量(i)会泄漏。

生成器

生成器

生成器是产生一列结果而不是单一值的函数。

David Beazley — A Curious Course on Coroutines and Concurrency

第三种创建迭代对象的方式是调用生成器函数。一个 生成器(generator) 是包含关键字yield的函数。值得注意，仅仅是这个关键字的出现完全改变了函数的本质：yield语句不必引发(invoke)，甚至不必可接触。但让函数变成了生成器。当一个函数被调用时，其中的指令被执行。而当一个生成器被调用时，执行在其中第一条指令之前停止。生成器的调用创建依附于迭代协议的生成器对象。就像常规函数一样，允许并发和递归调用。

当next被调用时，函数执行到第一个yield。每次遇到yield语句获得一个作为next返回的值，在yield语句执行后，函数的执行又被停止。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

>>> def f():
...   yield 1
...   yield 2
>>> f()                                   
<generator object f at 0x...>
>>> gen = f()
>>> gen.next()
1
>>> gen.next()
2
>>> gen.next()
Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
StopIteration

让我们遍历单个生成器函数调用的整个历程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

>>> def f():
...   print("-- start --")
...   yield 3
...   print("-- middle --")
...   yield 4
...   print("-- finished --")
>>> gen = f()
>>> next(gen)
-- start --
3
>>> next(gen)
-- middle --
4
>>> next(gen)                            
-- finished --
Traceback (most recent call last):
 ...
StopIteration

相比常规函数中执行f()立即让print执行，gen不执行任何函数体中语句就被赋值。只有当gen.next()被next调用，直到第一个yield部分的语句才被执行。第二个语句打印-- middle --并在遇到第二个yield时停止执行。第三个next打印-- finished --并且到函数末尾，因为没有yield，引发了异常。

当函数yield之后控制返回给调用者后发生了什么？每个生成器的状态被存储在生成器对象中。从这点看生成器函数，好像它是运行在单独的线程，但这仅仅是假象：执行是严格单线程的，但解释器保留和存储在下一个值请求之间的状态。^1

为何生成器有用？正如关于迭代器这部分强调的，生成器函数只是创建迭代对象的又一种方式。一切能被yield语句完成的东西也能被next方法完成。然而，使用函数让解释器魔力般地创建迭代器有优势。一个函数可以比需要next和__iter__方法的类定义短很多。更重要的是，相比不得不对迭代对象在连续next调用之间传递的实例(instance)属性来说，生成器的作者能更简单的理解局限在局部变量中的语句。

还有问题是为何迭代器有用？当一个迭代器用来驱动循环，循环变得简单。迭代器代码初始化状态，决定是否循环结束，并且找到下一个被提取到不同地方的值。这凸显了循环体——最值得关注的部分。除此之外，可以在其它地方重用迭代器代码。

双向通信

每个yield语句将一个值传递给调用者。这就是为何PEP 255引入生成器(在Python2.2中实现)。但是相反方向的通信也很有用。一个明显的方式是一些外部(extern)语句，或者全局变量或共享可变对象。通过将先前无聊的yield语句变成表达式，直接通信因PEP 342成为现实(在2.5中实现)。当生成器在yield语句之后恢复执行时，调用者可以对生成器对象调用一个方法，或者传递一个值 给 生成器，然后通过yield语句返回，或者通过一个不同的方法向生成器注入异常。

第一个新方法是send(value)，类似于next()，但是将value传递进作为yield表达式值的生成器中。事实上，g.next()和g.send(None)是等效的。

第二个新方法是throw(type, value=None, traceback=None)，等效于在yield语句处

1

raise type, value, traceback

不像raise(从执行点立即引发异常)，throw()首先恢复生成器，然后仅仅引发异常。选用单次throw就是因为它意味着把异常放到其它位置，并且在其它语言中与异常有关。

当生成器中的异常被引发时发生什么？它可以或者显式引发，当执行某些语句时可以通过throw()方法注入到yield语句中。任一情况中，异常都以标准方式传播：它可以被except和finally捕获，或者造成生成器的中止并传递给调用者。

因完整性缘故，值得提及生成器迭代器也有close()方法，该方法被用来让本可以提供更多值的生成器立即中止。它用生成器的__del__方法销毁保留生成器状态的对象。

让我们定义一个只打印出通过send和throw方法所传递东西的生成器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

>>> import itertools
>>> def g():
...     print '--start--'
...     for i in itertools.count():
...         print '--yielding %i--' % i
...         try:
...             ans = yield i
...         except GeneratorExit:
...             print '--closing--'
...             raise
...         except Exception as e:
...             print '--yield raised %r--' % e
...         else:
...             print '--yield returned %s--' % ans

>>> it = g()
>>> next(it)
--start--
--yielding 0--
0
>>> it.send(11)
--yield returned 11--
--yielding 1--
1
>>> it.throw(IndexError)
--yield raised IndexError()--
--yielding 2--
2
>>> it.close()
--closing--

注意： next还是__next__?

在Python 2.x中，接受下一个值的迭代器方法是next，它通过全局函数next显式调用，意即它应该调用__next__。就像全局函数iter调用__iter__。这种不一致在Python 3.x中被修复，it.next变成了it.__next__。对于其它生成器方法——send和throw情况更加复杂，因为它们不被解释器隐式调用。然而，有建议语法扩展让continue带一个将被传递给循环迭代器中send的参数。如果这个扩展被接受，可能gen.send会变成gen.__send__。最后一个生成器方法close显然被不正确的命名了，因为它已经被隐式调用。

链式生成器

注意： 这是PEP 380的预览(还未被实现，但已经被Python3.3接受)[^2]

比如说我们正写一个生成器，我们想要yield一个第二个生成器——一个子生成器(subgenerator)——生成的数。如果仅考虑产生(yield)的值，通过循环可以不费力的完成：

1
2
3

subgen = some_other_generator()
for v in subgen:
    yield v

然而，如果子生成器需要调用send()、throw()和close()和调用者适当交互的情况下，事情就复杂了。yield语句不得不通过类似于前一章节部分定义的try...except...finally结构来保证“调试”生成器函数。这种代码在PEP 380中提供，现在足够拿出将在Python 3.3中引入的新语法了：

1

yield from some_other_generator()

像上面的显式循环调用一样，重复从some_other_generator中产生值直到没有值可以产生，但是仍然向子生成器转发send、throw和close。

装饰器

总结

这个语言中令人激动的特性几乎充满歉意的，考虑到它可能没这么有用。

Bruce Eckel — An Introduction to Python Decorators

因为函数或类都是对象，它们也能被四处传递。它们又是可变对象，可以被更改。在函数或类对象创建后但绑定到名字前更改之的行为为装饰(decorator)。^4

“装饰器”后隐藏了两种意思——一是函数起了装饰作用，例如，执行真正的工作，另一个是依附于装饰器语法的表达式，例如，at符号和装饰函数的名称。

函数可以通过函数装饰器语法装饰：

1
2
3

@decorator             # ②
def function():        # ①
    pass

• 函数以标准方式定义。①
• 以@做为定义为装饰器函数前缀的表达式②。在 @ 后的部分必须是简单的表达式，通常只是函数或类的名字。这一部分先求值，在下面的定义的函数准备好后，装饰器被新定义的函数对象作为单个参数调用。装饰器返回的值附着到被装饰的函数名。

装饰器可以应用到函数和类上。对类语义很明晰——类定义被当作参数来调用装饰器，无论返回什么都赋给被装饰的名字。

在装饰器语法实现前(PEP 318)，通过将函数和类对象赋给临时变量然后显式调用装饰器然后将返回值赋给函数名，可以完成同样的事。这似乎要打更多的字，也确实装饰器函数名用了两次同时临时变量要用至少三次，很容易出错。以上实例相当于：

1
2
3

def function():                  # ①
    pass
function = decorator(function)   # ②

装饰器可以堆栈(stacked)——应用的顺序是从底到上或从里到外。就是说最初的函数被当作第一次参数器的参数，无论返回什么都被作为第二个装饰器的参数……无论最后一个装饰器返回什么都被依附到最初函数的名下。

装饰器语法因其可读性被选择。因为装饰器在函数头部前被指定，显然不是函数体的一部分，它只能对整个函数起作用。以@为前缀的表达式又让它明显到不容忽视(根据PEP叫在您脸上……：）)。当多个装饰器被应用时，每个放在不同的行非常易于阅读。

代替和调整原始对象

装饰器可以或者返回相同的函数或类对象或者返回完全不同的对象。第一种情况中，装饰器利用函数或类对象是可变的添加属性，例如向类添加文档字符串(docstring).装饰器甚至可以在不改变对象的情况下做有用的事，例如在全局注册表中注册装饰的类。在第二种情况中，简直无所不能：当什么不同的东西取代了被装饰的类或函数，新对象可以完全不同。然而这不是装饰器的目的：它们意在改变装饰对象而非做不可预料的事。因此当一个函数在装饰时被完全替代成不同的函数时，新函数通常在一些准备工作后调用原始函数。同样，当一个类被装饰成一个新类时，新类通常源于被装饰类。当装饰器的目的是“每次都”做什么，像记录每次对被装饰函数的调用，只有第二类装饰器可用。另一方面，如果第一类足够了，最好使用它因为更简单。^3

实现类和函数装饰器

对装饰器惟一的要求是它能够单参数调用。这意味着装饰器可以作为常规函数或带有__call__方法的类的实现，理论上，甚至lambda函数也行。

让我们比较函数和类方法。装饰器表达式(@后部分)可以只是名字。只有名字的方法很好(打字少，看起来整洁等)，但是只有当无需用参数定制装饰器时才可能。被写作函数的装饰器可以用以下两种方式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

>>> def simple_decorator(function):
...   print "doing decoration"
...   return function
>>> @simple_decorator
... def function():
...   print "inside function"
doing decoration
>>> function()
inside function

>>> def decorator_with_arguments(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       # in this inner function, arg is available too
...       print "doing decoration,", arg
...       return function
...   return _decorator
>>> @decorator_with_arguments("abc")
... def function():
...   print "inside function"
defining the decorator
doing decoration, abc
>>> function()
inside function

这两个装饰器属于返回被装饰函数的类别。如果它们想返回新的函数，需要额外的嵌套，最糟的情况下，需要三层嵌套。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

>>> def replacing_decorator_with_args(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       # in this inner function, arg is available too
...       print "doing decoration,", arg
...       def _wrapper(*args, **kwargs):
...           print "inside wrapper,", args, kwargs
...           return function(*args, **kwargs)
...       return _wrapper
...   return _decorator
>>> @replacing_decorator_with_args("abc")
... def function(*args, **kwargs):
...     print "inside function,", args, kwargs
...     return 14
defining the decorator
doing decoration, abc
>>> function(11, 12)
inside wrapper, (11, 12) {}
inside function, (11, 12) {}
14

_wrapper函数被定义为接受所有位置和关键字参数。通常我们不知道哪些参数被装饰函数会接受，所以wrapper将所有东西都创递给被装饰函数。一个不幸的结果就是显式参数很迷惑人。

相比定义为函数的装饰器，定义为类的复杂装饰器更简单。当对象被创建，__init__方法仅仅允许返回None，创建的对象类型不能更改。这意味着当装饰器被定义为类时，使用无参数的形式没什么意义：最终被装饰的对象只是装饰类的一个实例而已，被构建器(constructor)调用返回，并不非常有用。讨论在装饰表达式中给出参数的基于类的装饰器，__init__方法被用来构建装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

>>> class decorator_class(object):
...   def __init__(self, arg):
...       # this method is called in the decorator expression
...       print "in decorator init,", arg
...       self.arg = arg
...   def __call__(self, function):
...       # this method is called to do the job
...       print "in decorator call,", self.arg
...       return function
>>> deco_instance = decorator_class('foo')
in decorator init, foo
>>> @deco_instance
... def function(*args, **kwargs):
...   print "in function,", args, kwargs
in decorator call, foo
>>> function()
in function, () {}

相对于正常规则(PEP 8)由类写成的装饰器表现得更像函数，因此它们的名字以小写字母开始。

事实上，创建一个仅返回被装饰函数的新类没什么意义。对象应该有状态，这种装饰器在装饰器返回新对象时更有用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

>>> class replacing_decorator_class(object):
...   def __init__(self, arg):
...       # this method is called in the decorator expression
...       print "in decorator init,", arg
...       self.arg = arg
...   def __call__(self, function):
...       # this method is called to do the job
...       print "in decorator call,", self.arg
...       self.function = function
...       return self._wrapper
...   def _wrapper(self, *args, **kwargs):
...       print "in the wrapper,", args, kwargs
...       return self.function(*args, **kwargs)
>>> deco_instance = replacing_decorator_class('foo')
in decorator init, foo
>>> @deco_instance
... def function(*args, **kwargs):
...   print "in function,", args, kwargs
in decorator call, foo
>>> function(11, 12)
in the wrapper, (11, 12) {}
in function, (11, 12) {}

像这样的装饰器可以做任何事，因为它能改变被装饰函数对象和参数，调用被装饰函数或不调用，最后改变返回值。

复制原始函数的文档字符串和其它属性

当新函数被返回代替装饰前的函数时，不幸的是原函数的函数名，文档字符串和参数列表都丢失了。这些属性可以部分通过设置__doc__(文档字符串)，__module__和__name__(函数的全称)、__annotations__(Python 3中关于参数和返回值的额外信息)移植到新函数上，这些工作可通过functools.update_wrapper自动完成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>>> import functools
>>> def better_replacing_decorator_with_args(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       print "doing decoration,", arg
...       def _wrapper(*args, **kwargs):
...           print "inside wrapper,", args, kwargs
...           return function(*args, **kwargs)
...       return functools.update_wrapper(_wrapper, function)
...   return _decorator
>>> @better_replacing_decorator_with_args("abc")
... def function():
...     "extensive documentation"
...     print "inside function"
...     return 14
defining the decorator
doing decoration, abc
>>> function                           
<function function at 0x...>
>>> print function.__doc__
extensive documentation

一件重要的东西是从可迁移属性列表中所缺少的：参数列表。参数的默认值可以通过__defaults__、__kwdefaults__属性更改，但是不幸的是参数列表本身不能被设置为属性。这意味着help(function)将显式无用的参数列表，使使用者迷惑不已。一个解决此问题有效但是丑陋的方式是使用eval动态创建wrapper。可以使用外部external模块自动实现。它提供了对decorator装饰器的支持，该装饰器接受wrapper并将之转换成保留函数签名的装饰器。

综上，装饰器应该总是使用functools.update_wrapper或者其它方式赋值函数属性。

标准库中的示例

首先要提及的是标准库中有一些实用的装饰器，有三种装饰器：

• classmethod让一个方法变成“类方法”，即它能够无需创建实例调用。当一个常规方法被调用时，解释器插入实例对象作为第一个参数self。当类方法被调用时，类本身被给做第一个参数，一般叫cls。

  类方法也能通过类命名空间读取，所以它们不必污染模块命名空间。类方法可用来提供替代的构建器(constructor):

  1 2 3 4 5 6 7 8

  class Array(object): def __init__(self, data): self.data = data @classmethod def fromfile(cls, file): data = numpy.load(file) return cls(data)

这比用一大堆标记的__init__简单多了。

• staticmethod应用到方法上让它们“静态”，例如，本来一个常规函数，但通过类命名空间存取。这在函数仅在类中需要时有用(它的名字应该以_为前缀)，或者当我们想要用户以为方法连接到类时也有用——虽然对实现本身不必要。

• property是对getter和setter问题Python风格的答案。通过property装饰的方法变成在属性存取时自动调用的getter。

  1 2 3 4 5 6 7 8 9

  >>> class A(object): ... @property ... def a(self): ... "an important attribute" ... return "a value" >>> A.a <property object at 0x...> >>> A().a 'a value'

  例如A.a是只读属性，它已经有文档了：help(A)包含从getter方法获取的属性a的文档字符串。将a定义为property使它能够直接被计算，并且产生只读的副作用，因为没有定义任何setter。

  为了得到setter和getter，显然需要两个方法。从Python 2.6开始首选以下语法：

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

  class Rectangle(object): def __init__(self, edge): self.edge = edge @property def area(self): """Computed area. Setting this updates the edge length to the proper value. """ return self.edge**2 @area.setter def area(self, area): self.edge = area ** 0.5

通过`property`装饰器取代带一个属性(property)对象的getter方法，以上代码起作用。这个对象反过来有三个可用于装饰器的方法`getter`、`setter`和`deleter`。它们的作用就是设定属性对象的getter、setter和deleter(被存储为`fget`、`fset`和`fdel`属性(attributes))。当创建对象时，getter可以像上例一样设定。当定义setter时，我们已经在`area`中有property对象，可以通过`setter`方法向它添加setter，一切都在创建类时完成。

之后，当类实例创建后，property对象和特殊。当解释器执行属性存取、赋值或删除时，其执行被下放给property对象的方法。

为了让一切一清二楚[^5]，让我们定义一个“调试”例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
>>> class D(object):
...    @property
...    def a(self):
...      print "getting", 1
...      return 1
...    @a.setter
...    def a(self, value):
...      print "setting", value
...    @a.deleter
...    def a(self):
...      print "deleting"
>>> D.a                                    
<property object at 0x...>
>>> D.a.fget                               
<function a at 0x...>
>>> D.a.fset                               
<function a at 0x...>
>>> D.a.fdel                               
<function a at 0x...>
>>> d = D()               # ... varies, this is not the same `a` function
>>> d.a
getting 1
1
>>> d.a = 2
setting 2
>>> del d.a
deleting
>>> d.a
getting 1
1



属性(property)是对装饰器语法的一点扩展。使用装饰器的一大前提——命名不重复——被违反了，但是目前没什么更好的发明。为getter，setter和deleter方法使用相同的名字还是个好的风格。

一些其它更新的例子包括：

• functools.lru_cache记忆任意维持有限 参数：结果 对的缓存函数(Python 3.2)
• functools.total_ordering是一个基于单个比较方法而填充丢失的比较(ordering)方法(__lt__,__gt__，__le__等等)的类装饰器。

函数的废弃

比如说我们想在第一次调用我们不希望被调用的函数时在标准错误打印一个废弃函数警告。如果我们不像更改函数，我们可用装饰器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

class deprecated(object):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated()                    # doctest: +SKIP
    ... def f():
    ...     pass
    >>> f()                              # doctest: +SKIP
    f is deprecated
    """
    def __call__(self, func):
        self.func = func
        self.count = 0
        return self._wrapper
    def _wrapper(self, *args, **kwargs):
        self.count += 1
        if self.count == 1:
            print self.func.__name__, 'is deprecated'
        return self.func(*args, **kwargs)

也可以实现成函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

def deprecated(func):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated                      # doctest: +SKIP
    ... def f():
    ...     pass
    >>> f()                              # doctest: +SKIP
    f is deprecated
    """
    count = [0]
    def wrapper(*args, **kwargs):
        count[0] += 1
        if count[0] == 1:
            print func.__name__, 'is deprecated'
        return func(*args, **kwargs)
    return wrapper

while-loop移除装饰器

例如我们有个返回列表的函数，这个列表由循环创建。如果我们不知道需要多少对象，实现这个的标准方法如下：

1
2
3
4
5
6
7
8

def find_answers():
    answers = []
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
        answers.append(ans)
    return answers

只要循环体很紧凑，这很好。一旦事情变得更复杂，正如真实的代码中发生的那样，这就很难读懂了。我们可以通过yield语句简化它，但之后用户不得不显式调用嗯list(find_answers())。

我们可以创建一个为我们构建列表的装饰器：

1
2
3
4

def vectorized(generator_func):
    def wrapper(*args, **kwargs):
        return list(generator_func(*args, **kwargs))
    return functools.update_wrapper(wrapper, generator_func)

然后函数变成这样：

1
2
3
4
5
6
7

@vectorized
def find_answers():
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
        yield ans

插件注册系统

这是一个仅仅把它放进全局注册表中而不更改类的类装饰器，它属于返回被装饰对象的装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

class WordProcessor(object):
    PLUGINS = []
    def process(self, text):
        for plugin in self.PLUGINS:
            text = plugin().cleanup(text)
        return text

    @classmethod
    def plugin(cls, plugin):
        cls.PLUGINS.append(plugin)

@WordProcessor.plugin
class CleanMdashesExtension(object):
    def cleanup(self, text):
        return text.replace('—', u'\N{em dash}')

这里我们使用装饰器完成插件注册。我们通过一个名词调用装饰器而不是一个动词，因为我们用它来声明我们的类是WordProcessor的一个插件。plugin方法仅仅将类添加进插件列表。

关于插件自身说下：它用真正的Unicode中的破折号符号替代HTML中的破折号。它利用unicode literal notation通过它在unicode数据库中的名称(“EM DASH”)插入一个符号。如果直接插入Unicode符号，将不可能区分所插入的和源程序中的破折号。

更多例子和参考

• PEP 310(函数和方法装饰语法)
• PEP 3129(类装饰语法)
• http://wiki.python.org/moin/PythonDecoratorLibrary
• http://docs.python.org/dev/library/functools.html
• http://pypi.python.org/pypi/decorator
• Bruce Eckel
  • 装饰器I:介绍Python装饰器
  • Python装饰器II：装饰器参数
  • Python装饰器III：一个基于装饰器构建的系统

上下文管理器

上下文管理器是可以在with语句中使用，拥有__enter__和__exit__方法的对象。

1
2

with manager as var:
    do_something(var)

相当于以下情况的简化：

1
2
3
4
5

var = manager.__enter__()
try:
    do_something(var)
finally:
    manager.__exit__()

换言之，PEP 343中定义的上下文管理器协议允许将无聊的try...except...finally结构抽象到一个单独的类中，仅仅留下关注的do_something部分。

1.__enter__方法首先被调用。它可以返回赋给var的值。as部分是可选的：如果它不出现，enter的返回值简单地被忽略。
2.with语句下的代码被执行。就像try子句，它们或者成功执行到底，或者break，continue或return，或者可以抛出异常。无论哪种情况，该块结束后，__exit__方法被调用。如果抛出异常，异常信息被传递给__exit__，这将在下一章节讨论。通常情况下，异常可被忽略，就像在finally子句中一样，并且将在__exit__结束后重新抛出。

比如说我们想确认一个文件在完成写操作之后被立即关闭：

1
2
3
4
5
6
7
8
9

>>> class closing(object):
...   def __init__(self, obj):
...     self.obj = obj
...   def __enter__(self):
...     return self.obj
...   def __exit__(self, *args):
...     self.obj.close()
>>> with closing(open('/tmp/file', 'w')) as f:
...   f.write('the contents\n')